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ABSTRACT 

We describe a practical method for constructing axisymmetric two-integral galaxy 
models (with distribution functions of the form f{E,Lz), in which E is the or¬ 
bital energy, and is the vertical component of the angular momentum), based 
on Schwarzschild’s orbit superposition method. Other f{E, L2)-methods are mostly 
based on solving the Jeans equations or on finding the distribution function directly 
from the density, which often places restrictions on the shape of the galaxy. Here, 
no assumptions are made and any axisymmetric density distribution is possible. The 
observables are calculated (semi-)analytically, so that our method is faster than most 
previous, fully numerical implementations. Various aspects are tested extensively, the 
results of which apply directly to three-integral Schwarzschild methods. We show that 
a given distribution function can be reproduced with high accuracy and investigate 
the behaviour of the parameter that is used to measure the goodness-of-fit. Further¬ 
more, we show that the method correctly identihes the range of cusp elopes for which 
axisymmetric two-integral models with a central black hole do not exist. 

Key words: galaxies: elliptical and lenticular, cD - galaxies: kinematics and dynamics 
- galaxies: structure 


1 INTRODUCTION 

The fundamental quantity in galaxy dynamics is the dis¬ 
tribution function (DF), which specifies the distribution of 
the stars in a galaxy over position and velocity. By Jeans’ 
theorem, the DF is a function of the isolating integrals that 
are conserved by the corresponding potential (e.g., Lynden- 
Bell 1962; Binney 1982). Axisymmetric galaxies, which we 
will study here, conserve at least the two classical integrals 
of motion, the energy E and the vertical component of the 
angular momentum L^. The DF can generally not be mea¬ 
sured directly, but since observationally accessible quanti¬ 
ties, such as the projected density and the line-of-sight ve¬ 
locity, are simple moments of the DF, photometric and kine¬ 
matic observations can provide information on its proper¬ 
ties. In some cases, the part of a two-integral distribution 
function that is even in the velocities, f{E,L^), can be ob¬ 
tained analytically by using integral transforms to solve the 
relation between the DF and the intrinsic density p{R,z), 
where {R, ()>, z) are the usual cylindrical coordinates. The 
odd part can be found similarly from p{v,p}, where {v^} is 
the mean streaming velocity. These Laplace (Lynden-Bell 
1962; Lake 1981), Stieltjes (Hunter 1975b) and Laplace- 
Mellin (Dejonghe 1986) transforms have the drawback that 
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numerical implementation is difficult and that they require 
p and Vcj, to be explicit functions of the potential V and the 
cylindrical radius R. A more general formalism (Hunter & 
Qian 1993; Qian et al. 1995, hereafter Q95), the HQ-method, 
uses contour integration instead of integral transforms. This 
means it is simpler to implement, does not explicitly require 
p{R,V) and v,i,{R,V), but a suitable contour for the inte¬ 
gration has to be chosen, which may not always be at hand. 

Because of the drawbacks of these analytical for¬ 
malisms, numerical methods that can be applied to arbitrary 
potential-density pairs are more attractive. Various methods 
have been developed, many of which circumvent knowledge 
of the DF by solving the Jeans equations directly (van der 
Marel et al. 1994; Magorrian 1995), while others assume that 
the DF can be represented by a superposition of basis func¬ 
tions (Dehnen & Gerhard 1994; Kuijken 1995). Accurate 
estimates of the mass-to-light ratio M/L can be obtained 
in both ways (van der Marel 1991; Shaw et al. 1993; van 
den Bosch 1996). The predicted central black hole masses 
that are obtained with two-integral models (Magorrian et 
al. 1998) seem to overestimate the true values, but are still 
very useful to narrow down the parameter range for more 
sophisticated modeling (van der Marel et al. 1998, hereafter 
vdM98; Gebhardt et al. 2000; Bower et al. 2001). 

Numerical orbit integrations and observations of, e.g., 
the anisotropy of the stellar dispersions in the solar neigh¬ 
bourhood show that, in addition to E and a third inte- 
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gral is conserved for most orbits. In separable potentials, this 
third integral is exact and has a closed form (e.g., Kuzmin 
1956; de Zeeuw 1985), allowing a direct calculation of the DF 
(e.g., Dejonghe & de Zeeuw 1988). The most general fam¬ 
ily of these potentials corresponds to flattened mass mod¬ 
els with constant-density cores (de Zeeuw, Franx & Peletier 
1986), which provide a poor description of the inner regions 
of most galaxies, since these contain a central density cusp 
(Lauer et al. 1995). We conclude that many of the existing 
two-integral, as well as the analytical three-integral meth¬ 
ods, are applicable to a limited range of galaxy models. 

A flexible method for calculating numerical galaxy mod¬ 
els was designed by Schwarzschild (1979, 1982), who repre¬ 
sented the DF numerically by the occupation numbers in 
a superposition of building blocks. There are no restric¬ 
tions on the density or the potential, and no a priori as¬ 
sumptions have to be made about the shape or the degree 
of anisotropy of the galaxy. Schwarzschild’s original imple¬ 
mentation was aimed at reproducing a given triaxial density 
distribution and was subsequently applied to a large variety 
of galaxy models, from spherically and axially symmetric 
(Richstone & Tremaine 1984; Levison & Richstone 1985) to 
triaxial density distributions (e.g., Vietri 1986; Statler 1987; 
Schwarzschild 1993, Merritt & Fridman 1996; Siopis & Kan- 
drup 2000). More general versions of the method can also re¬ 
produce kinematical observations of spherically and axially 
symmetric galaxies that obey up to three integrals of mo¬ 
tion (Zhao 1996; Rix et al. 1997, hereafter R97; Cretton et al. 
1999, hereafter C99; Cretton et al. 2000; Hafner et al. 2000). 

In most implementations of Schwarzschild’s method, or¬ 
bits are used as individual building blocks. In non-separable 
potentials the third integral of motion is not known explic¬ 
itly, so that the orbital properties can only be obtained by 
solving the equations of motion numerically. This orbit in¬ 
tegration has to be carried out until the relevant quantities 
have averaged out and do not change upon a new time-step 
(Pfenniger 1984; Copin et al. 2000). Especially for orbits 
that are near-stochastic, or that have nearly commensurate 
fundamental frequencies, this condition is difficult to achieve 
and the orbit integration can be very time-consuming. Fur¬ 
thermore, an orbit fills a region in space that can have sharp 
edges, depending on the combination of integrals of mo¬ 
tion that it obeys. A superposition of a finite (and rela¬ 
tively small) number of orbits can therefore show artefacts 
that are caused by these edges. One way to solve this is 
to ‘blur’ orbits in phase-space randomly, another is to use 
building blocks that are smoother (e.g., Zhao 1996; Mer¬ 
ritt & Fridman 1996; R97; C99). Of particular interest are 
the so-called isotropic and two-integral components (ICs and 
TICs). Isotropic components, which are completely specified 
by their energy E, fill the volume inside the equipotential at 
E, while TICs, which are fixed once E and Lz are chosen, 
completely fill the corresponding zero-velocity curve (ZVC). 
The advantages of these building blocks is apparent: both 
the equipotential and the ZVC are smooth surfaces, and ICs 
and TICs can be considered as weighted combinations of all 
orbits with the same energy E (or E and Lz), including the 
irregular orbits. 

These building blocks can be used in Schwarzschild’s 
method in two ways. They can be added to every energy or 
every combination (E, Lz) of a (three-integral) orbit library, 
which will ‘automatically’ take care of the stochastic orbits 



R 

Figure 1. The meridional-plane cross section of the toroidal 
volume defined by the energy E and angular momentum Lz in 
an axisymmetric potential. Its boundary is the zero-velocity curve 
ZVC. When Lz = 0, the ZVC coincides with the equipotential 
(indicated by the drawn curve), and the hole around the z-axis 
disappears. When the angular momentum is maximal, the ZVC 
reduces to the circular orbit (indicated by the cross). 


(Zhao 1996; R97; C99). Alternatively, a fully isotropic or 
two-integral model can be constructed by using only ICs or 
TICs as building blocks. C99 suggested that this might be a 
practical way of constructing such models, and derived some 
analytic properties, but did not pursue these models further. 
We do so in this paper and concentrate on the two-integral 
case, since the properties of the isotropic components follow 
from those of the TICs by taking Lz = 0. 

In §^, we collect the properties of the TICS, and sum¬ 
marize how we use them in Schwarzschild’s method. The 
numerical aspects are discussed in detail in ^ We show in 
^ that our method can reproduce a model with a known 
analytical distribution function with high accuracy. We also 
study a set of mass models with a central black hole, and 
show that our method is able to detect when a self-consistent 
solution does not exist. We summarize our conclusions in 
Applications of our software to model the observed two- 
dimensional surface brightness and kinematics of nearby 
elliptical galaxies observed with SAURDN and STIS are de¬ 
scribed in two follow-up papers (Verolme et al. 2002, in 
preparation; McDermid et al. 2002, in preparation). 


2 TWO-INTEGRAL COMPONENTS 

We introduce the two-integral components (TICs) and cal¬ 
culate their observable properties, and then describe how we 
use them to construct two-integral models for galaxies. The 
properties of ICs follow by taking = 0. 
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2.1 Properties 

A TIC is completely specified by its energy Ej and angn- 
lar momentum Lzj (where j is a label for the TIC) and is 
defined by its distribution function 

^ f C5{E-Ej)5{Lz-Lz,j) inside ZVC, ,,, 

elsewhere, 

with C a constant that is used to normalize the mass of the 
TIC. The abbreviation ZVC denotes the zero-velocity curve 
in the {R, 2 :)-plane, defined as the locus of points for which 


i?kin = 14ft (f?, z) — E = e), 


( 2 ) 


where 

VMR,z) = V{R,z)-^, (3) 

is the effective potential. The ZVC defines a toroidal volume 
around the z-axis (Fig. with inner and outer radii Ri and 
i? 2 . The size of the central hole is determined by the value 
of the angular momentum: when Lz — 0 , there is no hole 
at all (i?i = 0 ) and the torus fills the equipotential at the 
given energy completely, so that the TIC reduces to the 
isotropic component (IC) at that energy. For the maximum 
value Lz = RcVc (with = RcdVss{Rc)/dRc the circular 
velocity at radius Rc) the torus reduces to the circular orbit 
at Ej, and Ri = R2 = Rc- 

The density pj {R, z) of a TIC is defined as the integral 
of the DF over all velocities, and is given by 


Pj{R,z) 


^ inside ZVC, 

R 

0 elsewhere. 


(4) 


The velocity moments are defined as the average of (powers 
of) the velocity components over all velocities. From sym¬ 
metry arguments and from the fact that / = f(E,Lz), the 
only nonzero first and second moments are 


(vl) 

( 4 ) 


7T C Lz,j 


(v^z) = ^[VMR,z) - Ej], 


(5) 


where these expressions are valid inside the ZVC. All mo¬ 
ments vanish outside the ZVC. 

We assume that the galaxy is observed at an inclination 
i, so that the projected coordinates are given by 


x' = y 

y' = —X cosi-|-« sini ( 6 ) 

z = X sin i + z cos i. 


where x' and y' are chosen in the plane of the sky, with 
y' along the projected minor axis of the galaxy, and z' is 
measured along the line of sight. 

The projected density Tij(x',y') is the integral of the 
intrinsic density along the line of sight 


^j{x',y')= J- 


ttC dz' 


-I- ( 2 ' sin i — y' cos i)^ 


(7) 


which can be evaluated analytically. The intrinsic density of 
a TIC is only defined for points that lie inside the toroidal 
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Figure 2. The line-of-sight crosses the ZVC in two (case b) or 
four (case a) points, depending on the location of (x',y'). The 
plot shows an edge-on configuration, so that the intersections are 
pairwize symmetric. 


volume defined by the ZVC, so that the integral (^ is defined 
on the z'-intervals delimited by the two or four points where 
the line-of-sight crosses the ZVC (Fig. 0). We give all expres¬ 
sions for four intersections, labeled as < 22 < 0 < 23 < 2 ^ 
The case of two intersections follows by taking 22 = 23 = 0. 
We obtain 


'^jix ,y ) — ^z^ — Sz' + Sz' “ ^z 

with 


Ej,' = - In 

sin i 


2 sin i(^R:,+DpJ 


where we have defined 

R'z,! = sj -I- ( 2 ' sin i — y' cos i)^ 

and 


( 8 ) 

fc = l,...,4, (9) 

( 10 ) 


= z sini — y cost. 


( 11 ) 


The limiting case of edge-on viewing (i = 0) is given by 


E.- = 


'xC z'l,. 


\/ x''^ + y'- 


fc = 


( 12 ) 


The velocity profile Cj{x', y', v^i), defined as the distri¬ 
bution of stars over velocity at the position {x', y'), is the 
integral of f{E, Lz) over the line of sight and the two veloc¬ 
ity components Vz,/ and Vyi in the plane of the sky. Since 
the distribution function of a TIC is the product of two <5- 
functions, the triple integral in the expression for jCj reduces 
to a single integral, 


^j{x',y',Vz') =C 


dz JejLzj , 


(13) 


where the remaining integration in z' is over the same inter¬ 
vals as in eq. (Q), and Je^l^ ^ is the Jacobian for the change 
of variables from (n^./, Vyi) to {Ej , Lz,j ). It can be written as 
(C99) 
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with 


1 

■y/^ + -B V^t + C v^, 


(14) 


A = 2 [V{x' ,y' ,z')-Ej] cos^ i+Dl,)-Llj, 

B = —2xLz,jSmi, (15) 

c = -r3, 

with R'^i and D^i from eqs (^^ and (§) , respectively. 

The integral of the velocity prohle over all allowed val¬ 
ues of v^i is equal to T,j{x' ,y'). The projected velocity mo¬ 
ments (n"/) (with n = 1 , 2 ,...) can be evaluated similarly, 
resulting in 


(«z') = 


fdz' 

1 ■ 

-C 

Ir^ 


= 

C 

fdz' 

\, 3B,, 

2 

J R'3 


3B= -4AC 

-, ^ T^/- 


with 

Uv^, = arcsin 


v^, R'^, + x'L^j sini/R'^, 
. V2 Bkin ^/x'^ cos^ i + Dl, 


(16) 


(17) 


and JejL^ j is given in eq. (^^. The values V- and v+ are 
the limits of the integration in v^i. 


2.2 Schwarzschild’s method 

The individual TICs can be used to model galaxies with 
Schwarzschild’s orbit superposition method. The procedure 
follows the usual steps (e.g., C99). First determine the den¬ 
sity distribution, either from a theoretical model or by de- 
projecting the observed projected surface density of a galaxy 
and choosing a mass-to-light ratio. Then calculate the cor¬ 
responding gravitational potential from Poisson’s equation. 
In this potential, dehne a collection of TICs by specifying 
the integrals {Ej,Lz,j). Jeans’ theorem guarantees that a 
sampling of the full range in energies and angular momenta 
results in a TIC library that is representative for this poten¬ 
tial. This can be achieved by sampling the radius on a grid 
of riB circular radii {Be}, chosen e.g. logarithmically on a 
range that includes > 99.9% of the mass. From this, the en¬ 
ergy grid {Ec} follows by calculating the energy of the circu¬ 
lar orbit at the {Be}. The angular momenta can then be 
chosen linearly between the minimum {Lz = 0 , orbits con¬ 
fined to the meridional plane) and maximum {Lz = ReVe, 
the circular orbit) at every Be. 

This choice of grids dehnes a library of nt = uetilz 
TICs. Schwarzschild’s method then determines the weighted 
sum of the TICs in this library that is closest to a data set 
Pi (generally consisting of both photometric and kinematic 
data). If Oij is the contribution of the j-th TIC to the i-th 
aperture, the problem reduces to solving for the weights 
in 

Nt 

'y ^ Ij Oij = Pi, i = 1 ,..., Np, (18) 

3 

where Np is the number of apertures for which observations 


are available. The 7 j, which determine the weight of each 
individual TIC in this superposition, are found by using a 
least-squares method, with the additional constraint that 
Ij > 0 (e.g., R97). Furthermore, 7 = '){E,Lz), since every 
TIC is labeled by a unique combination of E and Lz- The 
DF that corresponds to this set of weights can then be found 
from (Vandervoort 1984) 


DF(B,L.) = 


l{E,Lz) 

mE,Lz Ae,Lz 


(19) 


where mE,Lz is the mass of the TIC and Ae,Lz is the area 
of the cell around E and Lz in integral space. 

Since the matrix problem of eq. ( |l 8 | ) is often ill- 
conditioned, the weights and therefore the reconstructed 
DF may vary rapidly as a function of E and Lz- This 
can be avoided by adding regularisation constraints to 
the problem (e.g., Zhao 1996; R97). These constraints 
force the weights 'y{E,Lz) towards a smooth function 
of E and Lz by minimizing the n-th order derivatives 
9 " 7 (B, Lz)/dE'^, d"^{E, Lz)/dL^. The degree of smoothen- 
ing is determined by the order n and by the maximum value 
that the derivatives are allowed to have. 

The predicted observables that are obtained by calcu¬ 
lating the weighted sum of the TICs usually differ from the 
observed values Pi. There are various reasons for this: the 
DF of the galaxy may not be of the form f{E,Lz), model 
parameters such as inclination, stellar mass-to-light ratio or 
central black hole mass may not be chosen properly, the ob¬ 
servations can contain systematical errors, and the calcula¬ 
tions are done on hnite grids, which introduces discretization 
effects. To measure the discrepancy, we define a quality-of-fit 
parameter, , e-s 


2 


X 


E 


i = l 


( p:-P i 

\ APi 


2 


( 20 ) 


in which B* is the model prediction at the i-th aperture 
and APi is the (observational) uncertainty that is associated 
with Pi- 


3 NUMERICAL IMPLEMENTATION 

The numerical implementation of Schwarzschild’s method 
for three-integral axisymmetric models including kinematic 
constraints is described in detail in C99. Here we discuss 
the simpliheations for models containing TICs, which can 
be used to speed up the calculations. Readers who are more 
interested in applications, can skip to 


3.1 Storage grids 

All observables are calculated on grids that are adapted 
to the size of the galaxy and to the constraints that have 
to be reproduced. The meridional plane density is stored 
on a polar grid in the meridional plane, with a radial grid 
that is sampled logarithmically on the same range as the 
{Be} and a linear grid in polar angle between 0 and 7 r/ 2 . 
The projected density is calculated on a Cartesian grid in 
the sky plane, and the velocity profile is stored on a three- 
dimensional Cartesian grid in x',y' and u^/. The bin sizes 
of both grids are adapted to the resolution of the kinematic 
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Figure 3. A Z'VC{Ej ^ Lzj) in the meridional plane (the grey 
region), together with the cells of the meridional plane grid in 
radius and angle. The labels denote the method of integration 
(analytical or numerical, see text). 


data to be reproduced. The observables are averaged over 
these grids, so that a continuous observable Oj{x,v) is rep¬ 
resented by a discrete set of averaged grid entries 

(0.)= [ Ojix,v). (21) 

Jcell 

All observables are related to the DF through triple (in¬ 
trinsic density and velocity profile) or quadruple (projected 
density) integrals, two of which are trivial for TICs, because 
of their d-function behaviour. This means that the averages 
(Oj) can be written as no worse than quadruple integrals 
over the region enclosed by the cross-section of the ZVC 
and the grid cell. Depending on the relative position of the 
grid cell and the ZVC, the four integrals can be interchanged 
and the expressions can be simplified. In the following, we 
show how this is achieved for all observables. 


3.2 Meridional plane density 

The simple form of the meridional plane density (^) leads to 
considerable simplifications in the evaluation of the averaged 
meridional plane density (pj). The contribution of a general 
grid cell bordered by (r_,r+) and {9-,0+) is 

'■+ 

(pj) =2^C J dr J (22) 
r_ e_ 

For cells that are completely surrounded by the ZVC, which 
are relatively numerous for TICs with small energies (corre¬ 
sponding to large ZVCs), this integral can be evaluated and 
is given by 

{p^)=i.^C{rl-rl){e+-e.). (23) 

When, on the other hand, Ej is large, so that the ZVC is less 
extended, the cell is (much) larger than the ZVC. This means 

© 0000 RAS, MNRAS 000, sa 


Cell 

r_ e 

r+ e 

6 »_ e 

Eq. 

New range 

1 

[ 0 , Ri] 

[0,Ri] 

- 


0 


- 

2 

[ 0 , Ri] 

[Ri, R 2 ] 

- 





* 

[0,Ri] 

[R 2 , 00 } 

- 




[Ri,R 2 ] 

3 

[Ri, R 2 ] 

[Ri, R 2 ] 

[ 0 , 8zvc] 




['■-T-I-] 

4 

[Ri, R 2 ] 

[Ri, R 2 ] 

[dzvc, 




- 

5 

[Ri, R 2 ] 

[R 2 , 00 } 



21 


[r_, R 2 ] 


Table 1. For each cell in Fig. g (and for one cell that is not shown, 
labeled with *), the simplifications in the calculation of the aver¬ 
aged meridional plane density are given. The equation that has 
to be used is listed, or a value of 0 when the cell is located outside 
the ZVC. When eq. (§|) applies, the new integration limits are 
given. 


that the integration region in eq. ( p^ can be adjusted to 
match the ZVC size better, which speeds up the calculations 
considerably. Table ^ and Fig. show for which cells one can 
apply these two simplifications. Calculating all meridional 
plane densities for a library of 70 x 20 TICs takes a few 
minutes on a machine with a 1 GHz processor. 


3.3 Projected density 

The average of the projected density (^) over a cell in 
a Cartesian grid on the sky, bordered by {x-,x+) and 
{y-, V+), is given by (cf. eqs || and ^ 

(E,-) = (S,,) - (E.,) + (E.,) - (E.,), (24) 

with 


^ 
I’ smi 


+ “+ 

JdT Jdy' In 2 sini (/?(,/ +D^i^ , (25) 


(or, when i = 0, the integral over eq. (|lj)). Since z'f. = 
z'f,(x' ,y'), this expression cannot be simplified further with¬ 
out knowledge of the relation between z'l. and x' ,y' . Unfor¬ 
tunately, this relation is usually complicated and results in 
a very lengthy integrand, which rules out analytical evalua¬ 
tion. Therefore, we integrate eq. (|^ numerically by using 
a routine for more-dimensional numerical quadrature from 
the NAG-library, DOIFCF. These calculations are slightly 
more time-consuming than the ones for the meridional plane 
density; calculating the averaged projected densiti es fo r a li¬ 
brary of 1400 TICs on the machine mentioned in §3.2 takes 
of the order of thirty minutes. 


3.4 Velocity Profile 

The averaged Cm a. cell bordered by (x-,x+), {y-,y+) and 
(?;_,u+) equals 


y'+ "+ 2 '(")+ 

(jC) J y J J dz Jsj, I 


(26) 


z'(v)_ 


This integral is only defined on the region where the Ja¬ 
cobian exists, which is bordered by the points for which 
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Figure 4. The extreme velocities vi ^2 of the velocity profile as 
a function of the line-of-sight coordinate 2 ', at fixed and nonzero 
x' and y'. Only positive values of z' are drawn, a similar config¬ 
uration exists for z' < 0. The lower curve is ui, and the upper 
curve is V 2 - The horizontal dotted lines illustrate possible loca¬ 
tions of the velocity bin (u_,u+). The domain of integration in 
V and z' is the region enclosed by the curve and these bins (the 
grey area is an example; the new integration limits, 2 ^, 2 ^, 2 ', 2 ^ 
are indicated). 


A + Bv,,+Cvl, =0 ( 0 . We solve this equation for to 
find the following relation between the minimum and maxi¬ 
mum velocities of the velocity proHle, vi,V2, and z'■. 


-x' Lz,j sin i ± R'^, V 2 £’kin\/a:'^ cos^ i + D'^, 
viAz) = - ^72 -’ ( 27 ) 


with R'^i and D^i given in eqs (|l0[ ) and (|ll|), respectively. 
The solid curve in Fig. ^ shows these extreme velocities as 
a function of z', for arbitrary x' and y'. From eq. the 

ui,2 are directly proportional to Skin, which means they only 
exist inside the intervals [z^jZ^] and [23,24]. When 2' = 2', 
Skin = 0, and the cut-off velocities are 


Vl, 2 (z') = - 


x' L. 


z,j sini 


S'2 


( 28 ) 


where vi is the velocity at z' = Z2, 23 and V2 is reached when 
z! = 2:], 24. It follows that vi = —V2 at x' = 0 , so that £. is 
symmetric around u = 0. 

Since v is integrated over the interval v-,v+, the limits 
of integration in the inner integral of eq. (p6|), z'(v) -, z'(v)-1-, 
are determined by the cross-section of the velocity bin and 
the region enclosed by vi and V2 (the grey area in Fig. ^ 
is an example). The z'- integration cannot in general be 
carried out analytically. However, we can interchange the 
inner two integrals and, by using the quadratic dependence 
of the Jacobian on U2/, evaluate the integration over v^i: 




, 2 ^ ,H_ 


= C 


+ + 
fdx' jdy' Jdz 




( 29 ) 


Figure 5. A series of velocity profiles along the minor axis, for a 
double power law density distribution (js^ . The formalism that is 
outlined in the text is used to calculate £, but no average is made 
over the grid cells in x' and y' and very narrow bins in velocity 
are used to ensure that the resulting curve is smooth. 


where we have used eq. (|t^. This expression can be sim¬ 
plified even further, since for the values of z' where the bin 
extends outside the region enclosed by vi and V2 {v+ > V2, 
V- < vi), we can replace u_,+ with 111,2. In this case, we 
substitute eq. in eq. carry out the integration over 
z' and use eq. ( ^. The averaged velocity proHle then equals 

( 30 ) 

By calculating the intersections (if any) of V- and v+ with 
wi,2 (respectively, 25 < z'^ and z'^ < z'^), we separate the 
integration interval in z' in regions where we use eq. ( | 30 | ) 
and others where eq. (|^) applies. The averaged velocity 
profile in the grey area of Fig. ^ therefore consists of the 
averaged £ in the areas {23, 25} x {ui(2'), V2{z')} (area 3 a), 
Wb, z'a} X {u_, V2{z')} (area 3 b), {z'^, 4 } x ■!^+} (area 
3 c), and {z'^, z'^} x {ui(2'), u+} (area 3 d), which contribute, 
respectively, > ( 71 )^'.„2- i ^)and 

{£)^i^ ;,i Thfi total averaged velocity profile in this cell 

is then given by the sum of these four terms. A similar rea¬ 
soning applies to all other bins in Fig. ^ (Table ^). To check 
whether the resulting £ is smooth and resembles those cal¬ 
culated by C 99 , we have calculated a series of £’s along the 
minor axis (which means they are symmetric around w = 0). 
The result is shown in Fig. which can be compared to, e.g.. 
Fig. 5 of C 99 . 


3.5 PSF convolution 

In the previous subsections, we described how to obtain the 
average of the TIC observables over the cells of the storage 
grids. These quantities can only be compared with obser¬ 
vational data when point spread function (PSF) effects are 
taken into account correctly. This means that observation- 
ally accessible quantities such as the projected density and 

© 0000 RAS, MNRAS 000, 0-0 
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Cell 


n+ £ 

z',v^, € 

1 

(—oo, i)i] 

( —OO, Ui] 

{4.2^} X 

2 

(—oo, i)i] 

( —oo, Ui] 

{ 4 . 4 } X {r>i,r>-|-} 
{4-4} X {t'-.-y-i} 
{4-4} X {^^1,-y-i-} 

3 

(-oo,i)i] 

[0i, O 2 ] 

{4-4} X {'vi,v2} 
{4’4} X {v-,v2} 
{4-4} X {r>-,i^+} 
{4-4} X {^i--y-i-} 

* 

( —OO, Ul] 

[ 02 , 00 ) 

{4-4} X {r>i-r> 2 } 

4 

[■ 0 i,- 02 ] 

[ 01 , 02 ] 

{ 4 ’ 4 } X {v-,v2} 

{4-4} X {r>--i^+} 
{4-4} X {^i--!'-i-} 

5 

[■Oi, - 02 ] 

[ 02 , 00 ) 

{ 4 - 4 } X {'vi,v2} 

{4-4} X {v-,v2} 
{4-4} X {»’--^^+} 

{ 4 - 4 } X {'»i-'y-i-} 

6 

[f) 2 ,Oo) 

[ 02 , 00 ) 

{4-4} X {v-,v2} 
{4-4} X {i^--'«+} 

{ 4 - 4 } X {^^--^> 2 } 

7 

[- 02 , oo) 

[ 02 , 00 } 

{4-4} X {^>--^^ 2 } 


Table 2. The first column gives the label of the velocity bins in 
Fig. ^ (and for one bin that is not shown, denoted with *). The 
second and third column give the intervals in which the limits of 
the velocity bin, v— and are located. In the fourth column, 
the domains in {z' ,v^i) on which eq. ( ^9| ) and eq. are defined, 
are given. The total velocity profile is then obtained by adding 
the contributions of all these regions. 


the velocity profile have to be convolved with the appropri¬ 
ate PSF, which is usually described in terms of a (sum of) 
Gaussian functions. The extra two-dimensional integration 
that this convolution introduces can be carried out analyti¬ 
cally when it is combined with the average over the storage 
cell (see Appendix D of Q95). This means that an extra 
term, consisting of error functions, appears in the expres¬ 
sions for the observables that were given previously. 


4 EXAMPLES AND TESTS 

To test our method, we apply the numerical implementation 
to three theoretical potential-density pairs. We investigate 
whether our formalism rules out the same models as previous 
methods and determine the accuracy with which a theoreti¬ 
cal distribution function can be reproduced. This also allows 
us to study the behaviour of the (@) the regularisa- 
tion algorithm, which is similar to the one that is used in the 
three-integral Schwarzschild methods (C99). Applications to 
two-dimensional photometry and kinematics of galaxies are 
described in Verolme et al. and McDermid et al. 


4.1 Regularisation: Kuzmin-Kutuzov model 

We first test the ability of our method to reproduce a known 
f{E,Lz), given p and V, and investigate the effect of regu¬ 
larisation on the results. We use the density and potential 
of the Kuzmin-Kutuzov (1962) model 

M + c^) + c^) + + ^07 u 

^ ’ Itt {B? -\- z'^ + -\- 2uY ’ 

© 0000 RAS, MNRAS 000, E10 



Fig ure 6. The upper panel shows the Kuzmin-Kutuzov DF 
( |32[ ) as a function of energy E and fractional angular momen¬ 
tum Lz!Lz.mayi- The contour levels are chosen logarithmically 
and darker shading corresponds to larger values of the DF. The re¬ 
constructed DF ( [l9| ) for a two-integral Schwarzschild model that 
fits only the meridional plane masses is shown in the middle panel. 
There is a large discrepancy between the two functions. The lower 
panel shows the smooth DF that is found by including regulari¬ 
sation constraints in the fit. The contours are drawn at the same 
levels as in the top panel, showing that this DF reproduces the 
theoretical model very closely. The fractional difference between 
the two DFs on a range that includes > 99% of the mass is smaller 
than 10“^. 


V{R,z) 


_GM_ 

\/R^ + z'^ + -\-2T 


(31) 


with u — yj(? R^ + (z^ + c2). The density distribution is 

nearly spheroidal, and many of the properties of these mod¬ 
els, including the projected surface density, can be evaluated 
explicitly. 

The even part of the f{E, Lz) distribution function for 
the Kuzmin-Kutuzov models is given by (Dejonghe & de 
Zeeuw 1988; Batsleer & Dejonghe 1993) 


f{E,Lz) = 


^2^5/2 


VSa 


E 


( 32 ) 


Idt 


(1 - T) [(3 -h4x, - x?)(l - a;,)(l - T) -h 12t^] 

(1 - 2aEtVT^ +eLz tVThET 
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with 


Xe{t) = 2aEt 


1 + € yfz t 


(33) 


This integral can be expressed in terms of elementary func¬ 
tions, but this leads to an exceedingly lengthy expression, 
and it is best evaluated numerically. A grey-scale plot of the 
Kuzmin-Kutuzov DF as a function of energy and fractional 
angular momentum is shown in the upper panel of Fig. 6. 

We attempted to reproduce this DF by calculating a 
two-integral Schwarzschild model that hts the meridional 
plane mass (^). We used n_B = 70 and = 20, includ¬ 
ing Z/z = 0, and adopted 16 radial and 7 angular bins, like 
C99. With these parameters, the model density was fitted 
with very high accuracy (fractional difference smaller than 
0.001), but the TIC weights vary dramatically with energy 
and angular momentum, resulting in a very jagged DF (the 
middle panel in Fig. 6). Including regularisation constraints 
enforces the weights to vary smoothly with TIC index, re¬ 
sulting in a DF that reproduces the theoretical DF to within 
0.1%. This is illustrated in the lower panel of Fig. 6, which 
was obtained by requiring the same htting accuracy (AP in 
eq. for both the meridional plane masses and the regular¬ 
isation constraints. This means that both constraint types 
have equal influence on the Hnal model, so that the model 
DF is smooth, but still reproduces the density. 

In order to measure how fitting to the mass determines 
the kinematics of the model, we have calculated a theoretical 
velocity profile on the minor axis directly from the distribu¬ 
tion function eq. (^, and compared this with the predicted 
profile, given by the superposition of all the individual TIC 
velocity profiles. Fig. ^ demonstrates that the agreement 
between the two curves is significantly better than the typi¬ 
cal measurement error in observed velocity profiles (van der 
Marel & Franx 1993; Joseph et al. 2001). 


4.2 Models with a central black hole 

We now consider f{E, Lz) models with a central black hole, 
represented by a point mass with potential GM,/r. 

For spherically symmetric galaxies, there is a direct 
relation between the luminous density distribution and 
the isotropic distribution function f{E) (Eddington 1916). 
When the DF is required to be non-negative, this relation 
implies that no spherical, isotropic galaxies with a central 
black hole exist for cusp slopes shallower than —0.5 (e.g., 
Binney & Tremaine 1987). An indication that this restric¬ 
tion disappears when a higher degree of anisotropy is allowed 
is illustrated by the f{E,L^) (with L the angular momen¬ 
tum) that is obtained by placing all stars on circular orbits 
(Hunter 1975a). 

To investigate this effect for general anisotropic distri¬ 
bution functions f{E,Lz), Q95 applied the HQ-method to 
density distributions that are stratified on similar concen¬ 
tric spheroids. They considered models with a density prohle 
that follows a double power law, given by 





1 + -r 

\b ) 

\h 


where po and a are the central density and cusp slope, rrd = 
B? jq^ with q the intrinsic axis ratio, and b is the break 



Figure 7. The theoretical velocity profile for a point on the 
minor axis, calculated directly from eq. (^ (the solid line), and 
the superposition of all individual TIC velocity profiles (indicated 
with diamonds) for the Kuzmin-Kutuzov density distribution of 
eq. (^). The two curves agree to 0.1%. 

radius. Q95 found that all oblate {q < 1) double power law 
profiles correspond to a non-negative distribution function. 
However, adding a central black to these models limits the 
central cusp slope to a < —0.5, as in the spherical case. 
This restriction is independent of black hole mass, as long 
as M. > 0. 

De Bruijne et al. (1996) showed that this result is valid 
for a much larger family of cusped spheroidal density profiles 
with a central black hole; when the potential is spherically 
symmetric, the distribution function 

fiE,L\Lz)^Y.i-r^Y il^Y 9{E), (35) 

' ^ \ / \ -^max / 

corresponds to the density 

p(r,6») = ^D'"’^r-“sin2^6», (36) 

A 

in which the are given in terms of F-functions by eq. 

(17) of de Bruijne et al. The argument of a F-function must 
be positive, which implies that, near a black hole, where the 
potential falls off as 1/r (J = 1 in their notation), eqs disl) 
and (^^ exist only for 

a < A — p — i. (37) 

We see from eq. ( p5| ) that the limiting case of a two-integral 
DF f{E, Lz) is given by p = 0, so that eq. ( ^ reduces to 

a<A-i. (38) 

This means that density distributions that can be expanded 
in terms of sin6' must obey condition ( p^ on each of the 
terms of the expansion in order to have a physical DF near 
the black hole. The lowest-order term of this expansion, usu¬ 
ally with A = 0, places the most stringent condition on a, 
i.e. a < —0.5. 
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Figure 8. The quality-of-fit selfconsistent two-integral 

models as a function of cusp slope o;, at constant intrinsic flat¬ 
tening. The large for values in the interval [—0.5,0], and 

negligible for a < —0.5. This is caused by the non-existence of 
the self-consistent f(E^Lz) models with a > —0.5. 



A 


Figure 9. The quality of fit S’® ^ function of central cusp 
slope a and intrinsic flattening go i for a value of the flattening at 
infinity goo = 0.99. The dashed line shows the theoretical limit 
on the cusp slope for two-integral distribution functions and the 
crosses indicate combinations of a and go for which models were 
calculated. 


These results show that analytical methods are very 
powerful for deriving criteria for the existence of models. 
In contrast, Schwarzschild’s method always finds an orbit 
superposition that fits the constraints to some degree. The 
(120| ), determined from residnal errors, is the only mea¬ 
sure for the existence of the nnderlying distribntion fnnc- 
tion. A priori, it is not clear whether this parameter is able 
to reprodnce the sharp limits on the model parameters that 
are fonnd with analytical methods. For this reason, we con¬ 
structed two-integral models for the density profile ( |34| ) . We 
added a central point mass M, ~ 1% of the total galaxy 
mass, mimicking a central black hol e, a nd varied the cusp 
slope at constant flattening. As in §4.1, we took ue = 70 
and riL^ =20 (again including Lz = 0), and adopted 16 
radial and 7 angular bins. 

Fig. ^ shows the resulting function of a. for self- 

consistent two-integral models. While X < 10“® for values 
of a below —0.5, it increases rapidly as soon as ol increases 
above —0.5. This shows that our method recovers the models 
of Q95, and that our goodness-of-fit parameter X is a good 
discriminator of the (non)-existence of a given model. 

We also investigated the axisymmetric version of the 
triaxial models with central cusps that were introduced by 
de Zeeuw & Carollo (1996). The potential and density are 
given in Appendix A. The behaviour of these models is com¬ 
pletely determined by their axis ratios at small and large 
radii, go and goo, and the central cusp slope a. The density 
distribution is dimpled on the symmetry axis, and, for some 
choices of a, go, even becomes negative. In advance, we can 
therefore rule out a (small) range of parameters for which 
the DF does not exist. The physical range of go decreases 
when goo is decreased. At goo = 0.99, which is considered 
here, models with go < 0.05 show this behaviour. We study 
the range 0.1 < go < 1. 


© 0000 RAS, MNRAS 000, E10 


We added a central black hole with a mass ~ 1% of the 
total mass, and calculated self-consiste nt t wo-integral mod¬ 
els using the same grid settings as in §4.1, Fig. ^ shows the 
contours of constant X Si'S a function of a and go. As ex¬ 
pected, the X is zero for values a < —0.5 and rises rapidly 
above a = —0.5, indicating that self-consistent two-integral 
distribution functions do not exist at these values of the 
central cusp slope. However, for go < 0.4, the transition to 
> 0 occurs at values of a > —0.5, which means that 
these flatter models have an even larger freedom in possi¬ 
ble central cusp slopes. The points that contribute most to 
this nonzero X located in the center of the model, i.e., 
near the black hole, and experiments with a denser sampling 
in the integrals of motion give similar results. The limiting 
behavior of the density is (cf. eq ([A^) 


limp(r,6') = y{r) 

r —>0 


Xri + sin^ 9 


(39) 


so that the importance of the constant term in the expan¬ 
sion (^) decreases when go decreases. This may gradually 
relax the restriction on a in the numerical model. A similar 
effect has been found recently in oblate Sridhar-Touma mod¬ 
els (Jalali & de Zeeuw 2001), which resemble the power-law 
models in many respects. However, the values of go for which 
this occurs are much smaller than the observed flattenings 
of early-type galaxies. 


5 SUMMARY AND CONCLUSIONS 

We have presented an alternative to the existing analyti¬ 
cal and numerical methods for calculating two-integral dis¬ 
tribution functions. It is based upon Schwarzschild’s orbit 
superposition method and able to deal with arbitrary den- 
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sity profiles and galaxy potentials. Instead of orbital build¬ 
ing blocks, we use the so-called two-integral components, 
which are smoother and implicitly include stochastic orbits. 
Due to their delta-function behaviour, the observables can 
be calculated (semi-)analytically, which speeds up the cal¬ 
culations considerably. We checked that the method is able 
to reproduce a known combination of potential, density and 
distribution function and tested the regularisation algorithm 
in the process. This test shows that fitting the meridional 
plane masses alone is not enough to constrain the DF, while 
including a (modest) amount of regularisation is enough to 
smoothen the DF toward the theoretical curve. 

We have also tested the (|2^), which measures the 
quality-of-ht of the resulting model, against analytical in¬ 
vestigations. This demonstrates that a high value of can 
be taken as a sure indication that the model does not exist. 

These results show that our method is flexible and re¬ 
produces analytical results very closely. Additionally, many 
of the characteristics of the implementation that we used, 
such as the regularisation method and the x^-pa-rameter, 
are used in more general (three-integral) models. Regulari¬ 
sation is therefore also a very useful and necessary tool to 
smoothen the distribution function that is found by these 
methods, and the x^-parameter is a useful diagnostic to test 
the existence of the resulting models. Applications in which 
we calculate fully three-integral models (including TICs) of 
observed kinematics will follow in a subsequent paper. 
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APPENDIX A: A FAMILY OF 
AXISYMMETRIC MODELS 


The axisymmetric limit of the triaxial models introduced by 
de Zeeuw & Carollo (1996) is dehned by the potential 


V{r,9) = -u{r) -b i(3cos^6»- 

with 



u{r) = < 

. 

In 

r -b 1 

\( ^ "l 


2 -bg 

Vr-biy 

v{r) = 

AP+°‘ 

[r -b r2)4+“ ’ 


2+a 


- 1 


(Al) 


for a = —2 


for a 7 ^ —2 


(A2) 


where A and R 2 are constants, we have chosen GM = 1, and 
consider a > —3. The corresponding density distribution is 


9^,9) = g{r) 


f{r) 

9(r) 


- 1 -b I sin^ e 


with 


f{r) 

gp) 


(3 -b a)r°‘ 

47 r(r -b 1)4+“ 

Ar“ |^4r^ -b 4(6 -b a)r r 2 — a(5 -b a)r'^ 
47 r(r -b + 2 )®+“ 


(A3) 


(A4) 


The constants A and +2 are related to the axis ratios of 
the density distribution at small and large radii, qo and q^c, 
respectively, by 


A 


r2 


(3 -b g) (1 — qto) 

2{2 + qL) 

~(I ~ gto) (2 -b Qq °') Q (5 -b g) 

Ail-q^m^ + qL) 


(AS) 
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